1 COLORS plot

colors <- c("LMER" = "#DCA237", "Anova" = "#459B76")
lines <- c("Balanced" = "solid", "Unbalanced" = "dotted")

2 ESTIMATES VARIANCE REAL DATA

#############################
#PERFORMANCE PARAMETERS VALUE
#############################
data_PERF <- read.table(file=here::here("data", "DATACOMPLET_PERF.csv"), sep=",", header=TRUE)
data_PERF$Original_environment<-as.factor(data_PERF$Original_environment)
data_PERF <- data_PERF[data_PERF$Original_environment!="WT3",]
data_PERF <- data_PERF[data_PERF$Test_environment!="GF",]
data_PERF$Test_environment <- as.factor(data_PERF$Test_environment)
data_PERF$Population <- as.factor(data_PERF$Population)
data_PERF <- data_PERF[data_PERF$Test_environment!="Grape",] 
data_PERF <- droplevels(data_PERF)
data_PERF$IndicG0<-as.numeric(ifelse (as.character(data_PERF$Gen) == "G0", 1, 0))

nfruit_data <- nlevels(data_PERF$Original_environment)
nhab_data <- nlevels(data_PERF$Test_environment)
npop_per_fruit_data <- nlevels(data_PERF$Population)/nfruit_data
nrep_data <- mean(tapply(data_PERF$Nb_adults,list(data_PERF$Population,
                                                  data_PERF$Test_environment, 
                                                  data_PERF$Generation),length),na.rm = TRUE)
ntrial_data <- mean(data_PERF$Nb_eggs, na.rm = TRUE)
                              

#POISSON VARIANCE VALUE
mperf_poisson <- lme4::lmer(log(Nb_adults+1) ~ 1+ (1|Population) + 
                               (1|Original_environment:Test_environment) + 
                               (1|Original_environment:Test_environment:IndicG0),
                            data = data_PERF)

sdpop_datapoisson <- sqrt(unlist(lme4::VarCorr(mperf_poisson))[1])                               
sdfruithab_datapoisson <- sqrt(unlist(lme4::VarCorr(mperf_poisson))[3])
sdfruithab_ng_datapoisson <- sqrt(unlist(lme4::VarCorr(mperf_poisson))[2])   
sdfruithab_datapoisson <- sdfruithab_ng_datapoisson
sigma_datapoisson <- sigma(mperf_poisson)


#BINOMIAL VARIANCE VALUE
data_PERF_rate <- data_PERF[data_PERF$Nb_eggs>=data_PERF$Nb_adults,]
mperf_binomial <- lme4::lmer(asin(sqrt(Nb_adults/Nb_eggs)) ~ 1+ (1|Population) + 
                               (1|Original_environment:Test_environment) + 
                               (1|Original_environment:Test_environment:IndicG0),
                            data = data_PERF_rate)

sdpop_databinomial <- sqrt(unlist(lme4::VarCorr(mperf_binomial))[1])                               
sdfruithab_databinomial <- sqrt(unlist(lme4::VarCorr(mperf_binomial))[3])
sdfruithab_ng_databinomial <- sqrt(unlist(lme4::VarCorr(mperf_binomial))[2])   
sdfruithab_databinomial <- sdfruithab_ng_databinomial
sigma_databinomial <- sigma(mperf_binomial)




###########################
#PREFERENCE PARAMETERS VALUE
###########################
data_PREF <- read.table(file=here::here("data", "DATACOMPLET_PREF.csv"), sep=",", header=TRUE)
data_PREF$Original_environment<-as.factor(data_PREF$Original_environment)
data_PREF <- data_PREF[data_PREF$Original_environment!="WT3",]
dim(data_PREF)
## [1] 4428   12
data_PREF <- data_PREF[data_PREF$Test_environment=="Cranberry"|
                         data_PREF$Test_environment=="Cherry"|
                          data_PREF$Test_environment=="Strawberry",]
data_PREF$Test_environment <- as.factor(data_PREF$Test_environment)
data_PREF$Population <- as.factor(data_PREF$Population)
data_PREF$BoxID <- as.factor(data_PREF$BoxID)
data_PREF$Nb_eggs <- as.numeric(as.character(data_PREF$Nb_eggs))

data_PREF <- droplevels(data_PREF)
data_PREF$IndicG0<-as.numeric(ifelse (as.character(data_PREF$Gen) == "G0", 1, 0))

nfruit_data_box <- nlevels(data_PREF$Original_environment)
nhab_data_box <- nlevels(data_PREF$Test_environment)
npop_per_fruit_data_box <- nlevels(data_PREF$Population)/nfruit_data
nrep_data_box <- mean(tapply(data_PREF$Nb_eggs,list(data_PREF$Population,
                                                  data_PREF$Test_environment, 
                                                  data_PREF$Generation),length),na.rm = TRUE)


#POISSON BOX VARIANCE VALUE
mperf_poisson_box <- lme4::lmer(log(Nb_eggs+1) ~ 1+ (1|Population) + 
                               (1|Original_environment:Test_environment) + 
                               (1|Original_environment:Test_environment:IndicG0) +
                                  (1|BoxID),
                            data = data_PREF)

sdbox_datapoisson_box <- sqrt(unlist(lme4::VarCorr(mperf_poisson_box))[1])
sdpop_datapoisson_box <- sqrt(unlist(lme4::VarCorr(mperf_poisson_box))[2])                              
sdfruithab_datapoisson_box <- sqrt(unlist(lme4::VarCorr(mperf_poisson_box))[4])
sdfruithab_ng_datapoisson_box <- sqrt(unlist(lme4::VarCorr(mperf_poisson_box))[3])   
sdfruithab_datapoisson_box <- sdfruithab_ng_datapoisson_box
sigma_datapoisson_box <- sigma(mperf_poisson_box)

#NUMBER OF SIMULS

# Power Analyses
nb_simul <- 10000
#Save the same number of simuls for each SA value
nb_per_SA <- 1000
# 
# #Number of simuls for each SA non-genetic value
# tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)

#False positive rate 
nb_simul_fpr <- 100

3 BALANCED AND UNBALANCED DESIGN

3.1 Normal distribution

3.1.1 Balanced

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", design = "balanced",
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim1 <- as.data.frame(t(sim1))

sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", design = "balanced",
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim2 <- as.data.frame(t(sim2))


sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", design = "balanced",
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim3 <- as.data.frame(t(sim3))

# Concatenate value
sim <- rbind(sim1,sim2, sim3)

#Add indic SA
sim <- add_indic_sign(sim)

#Add number of simul for each SA value to keep only 100 simuls per SA value
sim <- add_sim_number_SA(sim)



## Local adaptation (genetic)
tapply(sim$IndicGen[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], length)
## -1.4 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2 
##    1    7   11   28   36   97  100  100  100  100  100  100  100  100  100  100 
##  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8 
##  100  100  100  100  100  100  100  100  100  100  100  100   53   24   12    1
tapply(sim$IndicGen_aov[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], mean)
##      -1.4      -1.2      -1.1        -1      -0.9      -0.8      -0.7      -0.6 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0800000 0.5700000 
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.5700000 0.5300000 0.6800000 0.6500000 0.6500000 0.6900000 0.6500000 0.7900000 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.8000000 0.8700000 0.8300000 0.8500000 0.8867925 0.8333333 0.8333333 1.0000000
## Local adaptation (non-genetic)
tapply(sim$IndicNonGen[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##      -1.5      -1.3      -1.2      -1.1        -1      -0.9      -0.8      -0.7 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -0.6      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.1800000 0.1600000 0.2300000 0.2700000 0.3700000 0.3200000 0.5500000 0.5200000 
##         1       1.1       1.2       1.3       1.4       1.5       1.6       1.7 
## 0.5100000 0.6700000 0.7600000 0.6700000 0.7600000 0.7076923 0.9032258 0.8750000 
##       1.8         2 
## 0.7500000 1.0000000
tapply(sim$IndicNonGen_aov[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##      -1.5      -1.3      -1.2      -1.1        -1      -0.9      -0.8      -0.7 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -0.6      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.2100000 0.1700000 0.2300000 0.2800000 0.3900000 0.3100000 0.5300000 0.5000000 
##         1       1.1       1.2       1.3       1.4       1.5       1.6       1.7 
## 0.5000000 0.6700000 0.7500000 0.6600000 0.7700000 0.7076923 0.9032258 0.8750000 
##       1.8         2 
## 0.7500000 1.0000000
tapply(sim$IndicGen, sim$SA_Gen, length)
## -1.4 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2 
##    1    7   11   28   36   97  142  260  399  623  815  980 1030 1100 1136 1102 
##  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8 
## 1516 2305 3063 3552 3458 2895 2129 1440  927  479  255  124   53   24   12    1
tapply(sim$IndicNonGen, sim$SA_NonGen, length)
## -1.5 -1.3 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1 
##    1    1    2    8   23   50   86  162  301  414  591  814  972 1083 1112 1054 
##  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.6  1.7 
## 1096 1471 2217 3183 3662 3440 2810 2173 1374  911  462  292  126   65   31    8 
##  1.8    2 
##    4    1
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
## -1.4 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2 
##    1    7   11   28   36   97  142  260  399  623  815  980 1030 1100 1136 1102 
##  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8 
## 1516 2305 3063 3552 3458 2895 2129 1440  927  479  255  124   53   24   12    1
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 0
data_prop_gen <- data_prop

#Plot
plot_genetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Normal data") +
   theme_LO_sober  + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_normal

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 0
data_prop_nongen <- data_prop

#Plot
plot_nongenetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_normal

3.1.2 Unbalanced (low)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              design = "unbalanced", disp_design = 1,
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              design = "unbalanced", disp_design = 1,
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              design = "unbalanced", disp_design = 1,
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim3 <- as.data.frame(t(sim3))

# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)




####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA

sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 1
data_prop_gen <- rbind(data_prop_gen,data_prop)

#Plot
plot_genetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Normal data") +
   theme_LO_sober  + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_normal

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 1
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

#Plot
plot_nongenetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_normal

3.1.3 Unbalanced (high)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              design = "unbalanced", disp_design = 4,
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              design = "unbalanced", disp_design = 4,
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              design = "unbalanced", disp_design = 4,
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)




####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA

sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 4
data_prop_gen <- rbind(data_prop_gen,data_prop)

#Plot
plot_genetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Normal data") +
   theme_LO_sober  + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_normal

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 4
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

#Plot
plot_nongenetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_normal

3.1.4 Plot

data_proportion_gen<-tidyr::gather(data_prop_gen, "method", "proportion", 2:3)
data_proportion_gen$method <- plyr::revalue(data_proportion_gen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_gen$design <- as.factor(data_proportion_gen$design)


#Plot genetic
plot_genetic_normal_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_gen, 
   aes(x = val_SA_genet, y = proportion, color = method, linetype = design)) + 
   geom_point(size = 2) + 
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Normal data") +
   ylim(0,1) +
   theme_LO_sober  + 
   ggtitle("Normal data") + 
   labs(color = "Methods", linetype = "Unbalanced\ndesign") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_genetic_normal_balancedvsunbalanced

#Data non-genetic 
data_proportion_nongen<-tidyr::gather(data_prop_nongen, "method", "proportion", 2:3)
data_proportion_nongen$method <- plyr::revalue(data_proportion_nongen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_nongen$design <- as.factor(data_proportion_nongen$design)

#Plot non-genetic
plot_nongenetic_normal_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_nongen, 
   aes(x = val_SA_nongenet, y = proportion, color = method, linetype = design)) + 
   geom_point(size = 2) + 
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   ggtitle("Normal data") +
   ylim(0,1) +
   theme_LO_sober  + 
   ggtitle("Normal data") + 
   labs(color = "Methods", linetype = "Unbalanced\ndesign") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_nongenetic_normal_balancedvsunbalanced

3.1.5 False positive rate analyses

3.1.5.1 No genetic and no non-genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "normal", design = "balanced",
              rho=0, rho_ng=0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
  
## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.02
mean(sim$IndicGen_aov)
## [1] 0.02
## Local adaptation (non-genetic)
mean(sim$IndicNonGen)
## [1] 0.01
mean(sim$IndicNonGen_aov)
## [1] 0.02

3.1.5.2 Genetic effects and no non-genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "normal", design = "balanced",
              rho=-0.1, rho_ng=0,  
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
tapply(sim$IndicGen, sim$SA_Gen, length)
## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 
##   3   4  14   6  18  17  13  12   7   4   2
tapply(sim$IndicGen_aov, sim$SA_Gen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.0000000 0.2500000 0.0000000 0.5000000 0.2222222 0.3529412 0.4615385 0.3333333 
##         1       1.1       1.2 
## 0.5714286 0.7500000 1.0000000
## Local adaptation (non-genetic)
mean(sim$IndicNonGen)
## [1] 0.02
mean(sim$IndicNonGen_aov)
## [1] 0.02
tapply(sim$IndicNonGen, sim$SA_Gen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.0000000 0.0000000 0.0000000 0.1666667 0.0000000 0.0000000 0.0000000 0.0000000 
##         1       1.1       1.2 
## 0.0000000 0.2500000 0.0000000
tapply(sim$IndicNonGen_aov, sim$SA_Gen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.0000000 0.0000000 0.0000000 0.1666667 0.0000000 0.0000000 0.0000000 0.0000000 
##         1       1.1       1.2 
## 0.1428571 0.0000000 0.0000000
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)
##      -0.7      -0.6      -0.5      -0.4      -0.3      -0.2      -0.1         0 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.6666667 0.0000000
tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)
##      -0.7      -0.6      -0.5      -0.4      -0.3      -0.2      -0.1         0 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7 
## 0.0000000 0.1111111 0.0000000 0.0000000 0.0000000 0.3333333 0.0000000

3.1.5.3 Non-genetic effects and no genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "normal",design = "balanced",
              rho=0, rho_ng=-0.1,  
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.1
mean(sim$IndicGen_aov)
## [1] 0.1
tapply(sim$IndicGen, sim$SA_Gen, mean)
##         -1       -0.9       -0.7       -0.6       -0.5       -0.4       -0.3 
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
##       -0.2       -0.1          0        0.1        0.2        0.3        0.4 
## 0.00000000 0.00000000 0.00000000 0.08333333 0.00000000 0.12500000 0.20000000 
##        0.5        0.6        0.7        0.9        1.2 
## 0.33333333 0.50000000 0.33333333 0.50000000 1.00000000
tapply(sim$IndicGen_aov, sim$SA_Gen, mean)
##         -1       -0.9       -0.7       -0.6       -0.5       -0.4       -0.3 
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
##       -0.2       -0.1          0        0.1        0.2        0.3        0.4 
## 0.00000000 0.00000000 0.00000000 0.08333333 0.00000000 0.12500000 0.20000000 
##        0.5        0.6        0.7        0.9        1.2 
## 0.33333333 0.50000000 0.33333333 0.50000000 1.00000000
tapply(sim$IndicGen, sim$SA_NonGen, mean)
##        0.2        0.3        0.4        0.5        0.6        0.7        0.8 
## 0.00000000 0.00000000 0.00000000 0.08333333 0.07692308 0.25000000 0.08333333 
##        0.9          1        1.1        1.2 
## 0.12500000 0.00000000 0.00000000 0.50000000
tapply(sim$IndicGen_aov, sim$SA_NonGen, mean)
##        0.2        0.3        0.4        0.5        0.6        0.7        0.8 
## 0.00000000 0.00000000 0.00000000 0.08333333 0.07692308 0.25000000 0.08333333 
##        0.9          1        1.1        1.2 
## 0.12500000 0.00000000 0.00000000 0.50000000
## Local adaptation (non-genetic)
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)
##        0.2        0.3        0.4        0.5        0.6        0.7        0.8 
## 0.00000000 0.33333333 0.08333333 0.08333333 0.42307692 0.37500000 0.58333333 
##        0.9          1        1.1        1.2 
## 0.50000000 0.25000000 0.66666667 1.00000000
tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)
##        0.2        0.3        0.4        0.5        0.6        0.7        0.8 
## 0.00000000 0.00000000 0.08333333 0.08333333 0.42307692 0.37500000 0.58333333 
##        0.9          1        1.1        1.2 
## 0.50000000 0.25000000 0.66666667 1.00000000

3.2 Poisson distribution

3.2.1 Balanced

rm(sim)
rm(data_proportion_gen, data_proportion_nongen, 
   sim_select_genetic, sim_select_nongenetic,
   data_prop, data_prop_gen, data_prop_nongen)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))

#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))

#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))

# Concatenate value
sim <- rbind(sim1,sim2, sim3)

#Add indic SA
sim <- add_indic_sign(sim)

#Add number of simul for each SA value to keep only 100 simuls per SA value
sim <- add_sim_number_SA(sim)


## Local adaptation (genetic)
tapply(sim$IndicGen[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], length)
## -2.2   -2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 
##    1    5    3    6   11   19   19   29   55   87   93  100  100  100  100  100 
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1 
##  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100 
##  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9    2  2.1  2.2  2.3  2.4  2.5  2.6 
##  100  100  100  100  100  100  100  100  100  100  100  100   67   25   19   17 
##  2.7  2.8  2.9    3  3.1 
##   11    2    5    2    1
tapply(sim$IndicGen_aov[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], mean)
##      -2.2        -2      -1.9      -1.8      -1.7      -1.6      -1.5      -1.4 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -1.3      -1.2      -1.1        -1      -0.9      -0.8      -0.7      -0.6 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0100000 0.0000000 0.0000000 0.0000000 0.0000000 0.0100000 
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0400000 0.0100000 0.0400000 0.1000000 0.0600000 0.1000000 0.1600000 0.1700000 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.2300000 0.3300000 0.5500000 0.6900000 0.8000000 0.8800000 0.9000000 0.8600000 
##       1.9         2       2.1       2.2       2.3       2.4       2.5       2.6 
## 0.9400000 0.9500000 0.9500000 0.9600000 0.9850746 1.0000000 1.0000000 1.0000000 
##       2.7       2.8       2.9         3       3.1 
## 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## Local adaptation (non-genetic)
tapply(sim$IndicNonGen[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##      -2.4      -2.1      -1.9      -1.8      -1.7      -1.6      -1.5      -1.4 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -1.3      -1.2      -1.1        -1      -0.9      -0.8      -0.7      -0.6 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0000000 0.0200000 0.0300000 0.0400000 0.0800000 0.0900000 0.1500000 0.2200000 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.2800000 0.2900000 0.4900000 0.5800000 0.6300000 0.6900000 0.7700000 0.8300000 
##       1.9         2       2.1       2.2       2.3       2.4       2.5       2.6 
## 0.8900000 0.9100000 0.9400000 0.9700000 0.9358974 0.9545455 0.8750000 1.0000000 
##       2.7       2.8       2.9         3 
## 1.0000000 1.0000000 1.0000000 1.0000000
tapply(sim$IndicNonGen_aov[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##      -2.4      -2.1      -1.9      -1.8      -1.7      -1.6      -1.5      -1.4 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -1.3      -1.2      -1.1        -1      -0.9      -0.8      -0.7      -0.6 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0100000 0.0100000 0.0100000 
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0000000 0.0400000 0.0200000 0.0400000 0.0800000 0.1100000 0.1600000 0.2100000 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.2900000 0.2900000 0.4700000 0.5700000 0.6400000 0.6900000 0.7800000 0.8300000 
##       1.9         2       2.1       2.2       2.3       2.4       2.5       2.6 
## 0.9000000 0.9000000 0.9400000 0.9700000 0.9358974 0.9318182 0.8750000 1.0000000 
##       2.7       2.8       2.9         3 
## 1.0000000 1.0000000 1.0000000 1.0000000
tapply(sim$IndicGen, sim$SA_Gen, length)
## -2.2   -2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 
##    1    5    3    6   11   19   19   29   55   87   93  151  199  237  332  439 
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1 
##  498  544  621  614  673  668  732  674  655  771  955 1257 1610 1831 2099 2184 
##  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9    2  2.1  2.2  2.3  2.4  2.5  2.6 
## 2194 2084 1793 1470 1198  951  675  530  373  254  136  121   67   25   19   17 
##  2.7  2.8  2.9    3  3.1 
##   11    2    5    2    1
tapply(sim$IndicNonGen, sim$SA_NonGen, length)
## -2.4 -2.1 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 
##    1    1    2    5    8   13   26   42   56   73  112  172  209  254  327  412 
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1 
##  506  554  617  622  699  684  685  642  639  739  905 1239 1622 1949 2116 2242 
##  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9    2  2.1  2.2  2.3  2.4  2.5  2.6 
## 2190 2009 1750 1486 1205  903  716  506  339  251  179  107   78   44   24   18 
##  2.7  2.8  2.9    3 
##   13    4    3    2
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
## -2.2   -2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 
##    1    5    3    6   11   19   19   29   55   87   93  151  199  237  332  439 
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1 
##  498  544  621  614  673  668  732  674  655  771  955 1257 1610 1831 2099 2184 
##  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9    2  2.1  2.2  2.3  2.4  2.5  2.6 
## 2194 2084 1793 1470 1198  951  675  530  373  254  136  121   67   25   19   17 
##  2.7  2.8  2.9    3  3.1 
##   11    2    5    2    1
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 0
data_prop_gen <- data_prop

#Plot
plot_genetic_poisson<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   ggtitle("Poisson data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_poisson

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 0
data_prop_nongen <- data_prop

#Plot
plot_nongenetic_poisson<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Poisson data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_poisson

3.2.2 Unbalanced (low)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
              design = "unbalanced", disp_design = 1,
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)


sim1 <- as.data.frame(t(sim1))

#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
              design = "unbalanced", disp_design = 1,
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)


sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
              design = "unbalanced", disp_design = 1,
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)


sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)




####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA

sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 1
data_prop_gen <- rbind(data_prop_gen,data_prop)

#Plot
plot_genetic_poisson_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_poisson_unbalanced

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 1
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

#Plot
plot_nongenetic_poisson_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Poisson data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_poisson_unbalanced

3.2.3 Unbalanced (high)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
              design = "unbalanced", disp_design = 4,
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
              design = "unbalanced", disp_design = 4,
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
              design = "unbalanced", disp_design = 4,
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)




####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA

sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 4
data_prop_gen <- rbind(data_prop_gen,data_prop)

#Plot
plot_genetic_poisson_unbalanced_high<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_poisson_unbalanced_high

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 4
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

#Plot
plot_nongenetic_poisson_unbalanced_high<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Poisson data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_poisson_unbalanced_high

3.2.4 Plot

data_proportion_gen<-tidyr::gather(data_prop_gen, "method", "proportion", 2:3)
data_proportion_gen$method <- plyr::revalue(data_proportion_gen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_gen$design <- as.factor(data_proportion_gen$design)


#Plot genetic
plot_genetic_poisson_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_gen, 
   aes(x = val_SA_genet, y = proportion, color = method, linetype = design)) + 
   geom_point(size = 2) + 
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   ylim(0,1) +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Unbalanced\ndesign") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_genetic_poisson_balancedvsunbalanced

#Data non-genetic 
data_proportion_nongen<-tidyr::gather(data_prop_nongen, "method", "proportion", 2:3)
data_proportion_nongen$method <- plyr::revalue(data_proportion_nongen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_nongen$design <- as.factor(data_proportion_nongen$design)

#Plot non-genetic
plot_nongenetic_poisson_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_nongen, 
   aes(x = val_SA_nongenet, y = proportion, color = method, linetype = design)) + 
   geom_point(size = 2) + 
   geom_line(size = 1) + 
   ylim(0,1) +
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   ggtitle("Poisson data") + 
   labs(color = "Methods", linetype = "Unbalanced\ndesign") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_nongenetic_poisson_balancedvsunbalanced

3.2.5 False positive rate analyses

3.2.5.1 No genetic and no non-genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho=0, rho_ng=0,   
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
  
## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.03
mean(sim$IndicGen_aov)
## [1] 0.03
## Local adaptation (non-genetic)
mean(sim$IndicNonGen)
## [1] 0.01
mean(sim$IndicNonGen_aov)
## [1] 0.01

3.2.5.2 Genetic effects and no non-genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho=-0.1, rho_ng=0,   
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
tapply(sim$IndicGen, sim$SA_Gen, mean)
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0000000 0.0000000 0.0000000 0.1250000 0.0000000 0.3333333 0.2500000 0.3076923 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.3636364 0.4444444 0.3750000 0.6666667 0.3750000 0.6000000 1.0000000 1.0000000 
##       1.9         2 
## 1.0000000 1.0000000
tapply(sim$IndicGen_aov, sim$SA_Gen, mean)
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0000000 0.0000000 0.0000000 0.1250000 0.0000000 0.3333333 0.2500000 0.3076923 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.3636364 0.4444444 0.3750000 0.6666667 0.3750000 0.6000000 1.0000000 1.0000000 
##       1.9         2 
## 1.0000000 1.0000000
## Local adaptation (non-genetic)
mean(sim$IndicNonGen)
## [1] 0.04
mean(sim$IndicNonGen_aov)
## [1] 0.04

3.2.5.3 Non-genetic effects and no genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "poisson",design = "balanced",
              rho=0, rho_ng=-0.1,  
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.12
mean(sim$IndicGen_aov)
## [1] 0.12
## Local adaptation (non-genetic)
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0000000 0.0000000 0.0000000 0.2500000 0.1250000 0.0000000 0.4666667 0.4615385 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.2000000 0.5000000 0.5000000 0.1666667 0.6666667 0.0000000 0.6666667 1.0000000 
##       1.9         2 
## 1.0000000 1.0000000
tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0000000 0.0000000 0.0000000 0.2500000 0.1250000 0.0000000 0.4000000 0.4615385 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.2000000 0.4000000 0.5000000 0.1666667 0.6666667 0.0000000 1.0000000 1.0000000 
##       1.9         2 
## 1.0000000 1.0000000

3.3 Binomial distribution

3.3.1 Balanced

rm(sim)
rm(data_proportion_gen, data_proportion_nongen, 
   sim_select_genetic, sim_select_nongenetic,
   data_prop, data_prop_gen, data_prop_nongen)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial", design = "balanced",
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial", design = "balanced",
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial", design = "balanced",
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim3 <- as.data.frame(t(sim3))


# Concatenate value
sim <- rbind(sim1,sim2, sim3)


#Add indic SA
sim <- add_indic_sign(sim)

#Add number of simul for each SA value to keep only 100 simuls per SA value
sim <- add_sim_number_SA(sim)


## Local adaptation (genetic)
tapply(sim$IndicGen[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], length)
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6 
##    1   21  100  100  100  100  100  100  100  100  100   18
tapply(sim$IndicGen_aov[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], mean)
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0300000 0.4500000 0.6100000 
##       0.3       0.4       0.5       0.6 
## 0.6900000 0.7600000 0.8100000 0.8333333
## Local adaptation (non-genetic)
tapply(sim$IndicNonGen[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0400000 0.2400000 
##       0.3       0.4       0.5       0.6       0.7 
## 0.3300000 0.6100000 0.7500000 0.7777778 1.0000000
tapply(sim$IndicNonGen_aov[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0100000 0.0500000 0.3000000 
##       0.3       0.4       0.5       0.6       0.7 
## 0.3300000 0.6200000 0.6700000 0.8333333 1.0000000
tapply(sim$IndicGen, sim$SA_Gen, length)
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6 
##    1   21  174  833 2397 3215 4678 9840 6730 1849  244   18
tapply(sim$IndicNonGen, sim$SA_NonGen, length)
##  -0.5  -0.4  -0.3  -0.2  -0.1     0   0.1   0.2   0.3   0.4   0.5   0.6   0.7 
##     1    13   172   891  2382  3186  4577 10066  6569  1851   273    18     1
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6 
##    1   21  174  833 2397 3215 4678 9840 6730 1849  244   18
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 0
data_prop_gen <- data_prop

#Plot
plot_genetic_binomial<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Binomial data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_binomial

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 0
data_prop_nongen <- data_prop

#Plot
plot_nongenetic_binomial<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Binomial data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_binomial

3.3.2 Unbalanced (low)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
              design = "unbalanced", disp_design = 1,
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
              design = "unbalanced", disp_design = 1,
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim2 <- as.data.frame(t(sim2))

#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
              design = "unbalanced", disp_design = 1,
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim3 <- as.data.frame(t(sim3))

# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)




####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA

sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 1
data_prop_gen <- rbind(data_prop_gen,data_prop)

#Plot
plot_genetic_binomial_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Binomial data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_binomial_unbalanced

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 1
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

#Plot
plot_nongenetic_binomial_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Binomial data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_binomial_unbalanced

3.3.3 Unbalanced (high)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
              design = "unbalanced", disp_design = 4,
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
              design = "unbalanced", disp_design = 4,
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
              design = "unbalanced", disp_design = 4,
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim3 <- as.data.frame(t(sim3))

# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)




####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA

sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 4
data_prop_gen <- rbind(data_prop_gen,data_prop)

#Plot
plot_genetic_binomial_unbalanced_high<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Binomial data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_binomial_unbalanced_high

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 4
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

#Plot
plot_nongenetic_binomial_unbalanced_high<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Binomial data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_binomial_unbalanced_high

3.3.4 Plot

data_proportion_gen<-tidyr::gather(data_prop_gen, "method", "proportion", 2:3)
data_proportion_gen$method <- plyr::revalue(data_proportion_gen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_gen$design <- as.factor(data_proportion_gen$design)


#Plot genetic
plot_genetic_binomial_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_gen, 
   aes(x = val_SA_genet, y = proportion, color = method, linetype = design)) + 
   geom_point(size = 2) + 
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Binomial data") +
   ylim(0,1) +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Unbalanced\ndesign") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_genetic_binomial_balancedvsunbalanced

#Data non-genetic 
data_proportion_nongen<-tidyr::gather(data_prop_nongen, "method", "proportion", 2:3)
data_proportion_nongen$method <- plyr::revalue(data_proportion_nongen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_nongen$design <- as.factor(data_proportion_nongen$design)

#Plot non-genetic
plot_nongenetic_binomial_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_nongen, 
   aes(x = val_SA_nongenet, y = proportion, color = method, linetype = design)) + 
   geom_point(size = 2) + 
   ylim(0,1) +
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   ggtitle("Binomial data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Unbalanced\ndesign") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_nongenetic_binomial_balancedvsunbalanced

3.3.5 False positive rate analyses

3.3.5.1 No genetic and no non-genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "binomial", design = "balanced",
              rho=0, rho_ng=0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
  
## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.02
mean(sim$IndicGen_aov)
## [1] 0.02
## Local adaptation (non-genetic)
mean(sim$IndicNonGen)
## [1] 0
mean(sim$IndicNonGen_aov)
## [1] 0.02

3.3.5.2 Genetic effects and no non-genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "binomial", design = "balanced",
              rho=-0.1, rho_ng=0,
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
tapply(sim$IndicGen, sim$SA_Gen, mean)
##       0.1       0.2       0.3       0.4 
## 0.0500000 0.2926829 0.4062500 0.5714286
tapply(sim$IndicGen_aov, sim$SA_Gen, mean)
##       0.1       0.2       0.3       0.4 
## 0.0500000 0.2926829 0.4062500 0.5714286
## Local adaptation (non-genetic)
mean(sim$IndicNonGen)
## [1] 0
mean(sim$IndicNonGen_aov)
## [1] 0

3.3.5.3 non-genetic effects and no genetic effects

### 
sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "binomial", design = "balanced",
              rho=0, rho_ng=-0.1,
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.07
mean(sim$IndicGen_aov)
## [1] 0.07
## Local adaptation (non-genetic)
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)
##       0.1       0.2       0.3       0.4 
## 0.0000000 0.2075472 0.3333333 0.6000000
tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)
##       0.1       0.2       0.3       0.4 
## 0.0000000 0.3018868 0.4074074 0.6000000

4 PLOT Balanced vs Unbalanced

POWER_ALL<-cowplot::plot_grid(plot_genetic_normal_balancedvsunbalanced,
                              plot_genetic_poisson_balancedvsunbalanced,
                              plot_genetic_binomial_balancedvsunbalanced,
                              plot_nongenetic_normal_balancedvsunbalanced,
                              plot_nongenetic_poisson_balancedvsunbalanced,
                              plot_nongenetic_binomial_balancedvsunbalanced,
                              ncol = 3, nrow = 2,
                              scale = c(1,  1))
POWER_ALL

#Save plot
name_plot<-paste0("Simul_powertest_ALLDATA_BALANCED_and_UNBALANCED",
      nb_simul,"nbsimul_",nb_per_SA,
      "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   POWER_ALL, base_height = 20/cm(1),
#                   base_width = 40/cm(1), dpi = 1200)

5 BOX EFFECT Absent vs Present

5.1 Non box

rm(sim)
rm(data_proportion_gen, data_proportion_nongen, 
   sim_select_genetic, sim_select_nongenetic,
   data_prop, data_prop_gen, data_prop_nongen)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))


# Concatenate value
sim <- rbind(sim1,sim2, sim3)

sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)

colors <- c("LMER" = "#DCA237", "Anova" = "#459B76")
lines_box <- c("Present" = "dotted", "Absent" = "solid")
####################################
###### Plot genetic
####################################
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$sd_box <- 0
data_prop_gen <- data_prop



################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$sd_box <- 0
data_prop_nongen <- data_prop

5.2 Box (low)

#With lower value of sdbox (close to sdpop): 
sdbox_datapoisson_box
##     BoxID 
## 0.4464253
sigma_datapoisson
## [1] 0.8727978
sdpop_datapoisson
## Population 
##   0.504277
sdfruithab_datapoisson
## Original_environment:Test_environment:IndicG0 
##                                     0.8111606
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 0.5,
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 0.5,
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = sdbox_datapoisson_box,
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)

sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)


####################################
###### Plot genetic
####################################
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$sd_box <- 0.5
data_prop_gen <- rbind(data_prop_gen,data_prop)


################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$sd_box <- 0.5
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

5.3 Box (high)

#With lower value of sdbox (close to sdpop): 
sdbox_datapoisson_box
##     BoxID 
## 0.4464253
sigma_datapoisson
## [1] 0.8727978
sdpop_datapoisson
## Population 
##   0.504277
sdfruithab_datapoisson
## Original_environment:Test_environment:IndicG0 
##                                     0.8111606
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 1.5,
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 1.5,
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 1.5,
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)

sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)


####################################
###### Plot genetic
####################################
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))

data_prop$sd_box <- 1.5
data_prop_gen <- rbind(data_prop_gen,data_prop)



################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$sd_box <- 1.5
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

5.4 Plot

#################################
### Common plot
#################################
data_proportion_gen<-tidyr::gather(data_prop_gen, "method", "proportion", 2:3)
data_proportion_gen$method <- plyr::revalue(data_proportion_gen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_gen$sd_box <- as.factor(data_proportion_gen$sd_box)


#Plot genetic
plot_genetic_poisson_ALL_box<-ggplot2::ggplot(data = data_proportion_gen, 
   aes(x = val_SA_genet, y = proportion, color = method, linetype = sd_box)) + 
   geom_point(size = 2) + 
   ylim(0,1) +
   geom_line(size = 1, alpha = 0.7) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Box effect\n(sd box)") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_genetic_poisson_ALL_box

#Data non-genetic 
data_proportion_nongen<-tidyr::gather(data_prop_nongen, "method", "proportion", 2:3)
data_proportion_nongen$method <- plyr::revalue(data_proportion_nongen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_nongen$sd_box <- as.factor(data_proportion_nongen$sd_box)



#Plot non-genetic
plot_nongenetic_poisson_ALL_box<-ggplot2::ggplot(data = data_proportion_nongen, 
   aes(x = val_SA_nongenet, y = proportion, color = method, linetype = sd_box)) + 
   geom_point(size = 2) + 
   ylim(0,1) +
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Box effect\n(sd box)") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_nongenetic_poisson_ALL_box

power_poisson_box<-cowplot::plot_grid(plot_genetic_poisson_ALL_box,
                                      plot_nongenetic_poisson_ALL_box,
                          labels=c("A", "B"), 
                          label_size = 15,
                          ncol =1, nrow = 2,
                     #    hjust = 0, vjust = 1,
                          scale = c(1,  1))
power_poisson_box

#Save plot
name_plot<-paste0("Simul_powertest_poisson_boxeffect_",
      nb_simul,"nbsimul_",nb_per_SA,
      "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   power_poisson_box, base_height = 20/cm(1), base_width = 14/cm(1), dpi = 1200)

6 SIMUL WITH REAL UNBALANCED DESIGN

6.1 Dataset

#unbalanced dataset
data_PERF <- read.table(file=here::here("data", "DATACOMPLET_PERF.csv"), sep=",", header=TRUE)
data_PERF$Original_environment<-as.factor(data_PERF$Original_environment)
data_PERF <- data_PERF[data_PERF$Original_environment!="WT3",]
data_PERF <- data_PERF[data_PERF$Test_environment!="GF",]
data_PERF <- data_PERF[data_PERF$Test_environment!="Grape",] #mvrnom doesn't work with unbalanced matrix
data_PERF <- droplevels(data_PERF)
tapply(data_PERF$Population, list(data_PERF$Generation, data_PERF$Population), length)
##    Blackberry31 Blackberry32 Blackberry33 Blackberry34 Blackberry35
## G0           38           42           25            4           NA
## G2           27           39           44           11           46
##    Blackberry36 Blackberry37 Blackberry38 Blackberry39 Blackberry40
## G0            8           15           21           10           15
## G2           17           33            5           24           30
##    Blackberry43 Blackberry44 Blackberry45 Cherry103 Cherry104 Cherry3 Cherry47
## G0            8           12           43         3        39       9       36
## G2           10           30            4        19        16       9       39
##    Cherry50 Cherry51 Cherry52 Cherry6 Cherry7 Strawberry42 Strawberry44
## G0       36       48       45      36       6           15           21
## G2       11        1        4       7      37           46           30
##    Strawberry53
## G0            6
## G2           63

6.2 Normal distribution

####################################
###### Simulation
####################################
#Simul data 
sim <-  sapply(1:nb_simul, simul_fitnessdata_unbalanced, distrib = "normal", 
               unbalanced_dataset = data_PERF, 
               sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, 
              rho = -0.05, rho_ng = -0.05, 
              sigma = 0.5)
             

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)

####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5 
##   10  107  496 1127 1633 1919 1672 1307  868  464  234  106   37   19    1
#Save the same number of simuls for each SA value
nb_per_SA = 50
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))


#Plot
plot_genetic_normal_unbalanced_real<-ggplot2::ggplot(data = data_prop, 
                                                     aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Normal data (unbalanced experimental design)") +
   theme_LO_sober +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_normal_unbalanced_real

################################# 
### Plot non-genetic
#################################
#Number of simuls for each SA non-genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.7 
##    4  111  493 1054 1759 1912 1695 1232  853  451  267   98   52   13    5    1
#Save the same number of simuls for each SA value
nb_per_SA 
## [1] 50
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))


#Plot
plot_nongenetic_normal_unbalanced_real<-ggplot2::ggplot(data = data_prop, 
                                                   aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Normal data (unbalanced experimental design)") +
   labs(color = "Methods") +
   scale_color_manual(values = colors)
plot_nongenetic_normal_unbalanced_real

6.3 Poisson distribution

####################################
###### Simulation
####################################
#Simul data 
#nb_simul <- 1000
sim <- sapply(1:nb_simul, simul_fitnessdata_unbalanced, distrib = "poisson", 
              unbalanced_dataset = data_PERF, sdpop = 0.5, 
              sdfruithab = 1, sdfruithab_ng = 1, rho = -0.05, rho_ng = -0.05, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)

####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9   2 2.1 
##   3  10  50 121 253 404 558 709 819 939 955 915 814 777 652 531 465 313 221 171 
## 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.1 
## 119  63  60  36  15  12  11   3   1
#Save the same number of simuls for each SA value
# nb_per_SA=10
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))


#Plot
plot_genetic_poisson_unbalanced_real<-ggplot2::ggplot(data = data_prop, 
                                                      aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data (unbalanced experimental design)") +
   theme_LO_sober  + 
   labs(color = "Methods") +
   scale_color_manual(values = colors)
plot_genetic_poisson_unbalanced_real

################################# 
### Plot non-genetic
#################################
#Number of simuls for each SA non-genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9   2 2.1 
##   2  11  49 129 240 379 522 727 885 969 951 905 845 743 663 496 405 325 224 183 
## 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9   3 3.3 
## 127  78  53  35  26  15   5   3   4   1
#Save the same number of simuls for each SA value
nb_per_SA 
## [1] 50
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))


colors<-c("LMER" = "#DCA237", "Anova" = "#459B76")
#Plot
plot_nongenetic_poisson_unbalanced_real<-ggplot2::ggplot(data = data_prop, 
                                                   aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Poisson data (unbalanced experimental design)") +
   labs(color = "Methods") +
   scale_color_manual(values = colors)
plot_nongenetic_poisson_unbalanced_real

6.4 Binomial distribution

####################################
###### Simulation
####################################
#Simul data 
#nb_simul <- 1000
sim <- sapply(1:nb_simul, simul_fitnessdata_unbalanced, distrib = "binomial", 
              unbalanced_dataset = data_PERF, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, 
              rho = -0.05, rho_ng = -0.05, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)

####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5 
##   10  107  496 1127 1633 1919 1672 1307  868  464  234  106   37   19    1
#Save the same number of simuls for each SA value
# nb_per_SA=10
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))


#Plot
plot_genetic_Binomial_unbalanced_real<-ggplot2::ggplot(data = data_prop, 
                                                  aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Binomial data (unbalanced experimental design)") +
   theme_LO_sober  + 
   labs(color = "Methods") +
   scale_color_manual(values = colors) 
plot_genetic_Binomial_unbalanced_real

################################# 
### Plot non-genetic
#################################
#Number of simuls for each SA non-genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.7 
##    4  111  493 1054 1759 1912 1695 1232  853  451  267   98   52   13    5    1
#Save the same number of simuls for each SA value
nb_per_SA 
## [1] 50
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))


colors<-c("LMER" = "#DCA237", "Anova" = "#459B76")
#Plot
plot_nongenetic_Binomial_unbalanced_real<-ggplot2::ggplot(data = data_prop, 
                                                   aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Binomial data (unbalanced experimental design)") +
   labs(color = "Methods") +
   scale_color_manual(values = colors)
plot_nongenetic_Binomial_unbalanced_real

6.5 All plots unbalanced

power_plot_all_unbalanced_real<-cowplot::plot_grid(plot_nongenetic_normal_unbalanced_real,
                                   plot_nongenetic_poisson_unbalanced_real,
                                   plot_nongenetic_Binomial_unbalanced_real,
                                   plot_genetic_normal_unbalanced_real,
                                   plot_genetic_poisson_unbalanced_real,
                                   plot_genetic_Binomial_unbalanced_real,
                          #labels=c("A", "B"), 
                          #label_size = 15,
                          ncol =3, nrow = 2,
                     #    hjust = 0, vjust = 1,
                          scale = c(1,  1))
power_plot_all_unbalanced_real

#Save plot
name_plot<-paste0("Simul_powertest_ALLDATA_UNBALANCED_real",
      nb_simul,"nbsimul_",nb_per_SA,
      "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   power_plot_all_unbalanced_real, base_height = 20/cm(1),
#                   base_width = 42/cm(1), dpi = 1200)
#